function [L,U]=mylu(A)
[n,n]=size(A);
for i=1:n   
    L(i,i)=1;
end
U(1,1)=A(1,1)/L(1,1);
if L(1,1)*U(1,1)==0
	fprintf('Factorization impossible')
else
	for j=2:n
		U(1,j)=A(1,j)/L(1,1);
		L(j,1)=A(j,1)/U(1,1);
	end
end
for i=2:n-1
	sum=0;
	for k=1:i-1
		sum=sum+L(i,k)*U(k,i);
	end
	U(i,i)=(A(i,i)-sum)/L(i,i);
	if L(i,i)*U(i,i)==0
		fprintf('Factorization impossible')
	else
		for j=i+1:n
			h=0;
			s=0;
			for k=1:i-1
				h=h+L(i,k)*U(k,j);
				s=s+L(j,k)*U(k,i);
			end
			U(i,j)=1/L(i,i)*(A(i,j)-h);
			L(j,i)=1/U(i,i)*(A(j,i)-s);
		end
	end
end
sum=0;
for k=1:n-1
	sum=sum+L(n,k)*U(k,n);
	U(n,n)=(A(n,n)-sum)/L(n,n);
end
